TI89 program “DIY” to have position with sun

Relative to the stars, it's still more interesting to calculate position with the sun; during the day the most of time the skyline is clear and easier findable, making sextant's use easier and more accurate, with a wider range of practices's time . As the program does with the stars, we will first search to “transform” all data of our time measurement (hours, minutes, seconds, month, years) in one lonely number. this last one will help us to unroll calculations reminder which is no more than functions woking with time. For this beginning we could have used the same as the stars program, but just for fun, This time we will take account of the Gregorian calendar, and get rid of julian calendar!

1) First part: how to calculate time by using gregorian calendar

Recall: in julian calendar the average year lasts 365.25 days, so a leap year comes each 4 years (2008, 2012, 2016, 2020…). on the calendar it means that we add up a day at the end of february each 4 years (29 days instead of 28).

Summary, on the way to have a calendar year period close to reality, in fact a tropic year period (nearly 365,2422 days), the calculations method has been improved to reduce slowly in time the julian calendar size. With the gregorian calendar an average year lasts 365,2425 days. to reach that point, we keep the julian's calendar features , but each secular year we consider it as a bissextile year only if it could be divided by 400. for our programmation we will then consider a year as bissextile if it is 4 years after an another bissextile's year. Then for secular's years which could be divided by 100 we select and we will accept them as bissextile if they are dividable too by 400. for exemple according to gregorian's calendar 2000 and 2400 are bissextile, but 2500 and 2600 are not bissextile.

here you will find the beginning:

 

 

we use the loop “while” and "int" function(). Loop “while” asks to the program to calculate something until it reaches one or several conditions, and function “int()” ne garde que la partie entière d’un nombre, for exemple int(2.7)=2. The number at the end of the page 2457388.5 is the julian's day 01/01/2016 at 00:00.

2) second part: how to calculate with time declination and hour angle 's sun.

As we have seen in “programming a TI89 for navigation”, there are many ways to calculate those parameters, more or less accurate on a short or long period, with more or less difficulties…to have accuracy in results on a very long period, it needs a lot of parameters to consider, and in addition planets have a chaotic behavior on wide periods of time…we will then go ahead all this and keep just all we need for hour needs, a program usefull and relevant for a life period and enough accurate to make a point in navigation and not too bothering to write!

a) hour angle:

if we consider things on the easiest side, we follow the fact that earth is spinning in 24 hours, giving illusion of sun's motion around earth moving at 360°/24H=15°per hour, with sun's passage at greenwich's meridian at 12H UTC. By keeping this writting AHG=15*(l*24) would be enough to calculate greenwich hour angle(AHG), with a greenwich's meridian noon as origin data for time. to find our local hour angle (AHL), only substract our dead longitude from AHG would be enough(G).

however it's inadequate, with this simple method at the end we find high gaps more and more 3 wide regards to the truth! in fact apparent sun's motion regards to earth has not a steadily speed because of the elliptic shape of the earth's way (remember the first and second Kepler's laws), and because of the rotation's earth axis inclined too. Fortunately, Gaps reiterates almost the same years by years, allowing quite easily to calculate good results each time. For this we will work with the time's equation. as our friend wikibobby says this equation matches for a precised time difference between average solar time and true solar time. the average solar time is calculation written in the paragraph just above (sun is moving with a constant velocity) and the true solar time will be calculated with those equations given by wikibobby (take care of earth elliptic's way and earth obliquity). Now in adequation with our program's beginning it gives us this:

The mod function() here is used to get a result always being between 0 and 360 degrees (This function provides the euclidian division's remainder depending on dividend and divider choosen), otherwise we would end up with huge angles without meaning. We could also use the loop “for” but for that the calculator would take more time to calculate. There is 180 after mod() because the program begins at midnight 1/1/2016 (see prog's first part), so this time the average hour angle of the sun with greenwich is 180 degrees (if the prog would have begun at noon, we wouldn't have added up 180 but nothing instead).

For the equations of time, coefficients used are easily noteworthy:

357.5291= average anomaly 01/01/2000 to 12h.

0.98560028= it's the angular's velocity of the average's rotation of earth around the sun between the passage of 2 perihelias, so it's 360 divided by the anomalistic year.

2451545= it's the julian's day of 01/01/2000 to 12h.

1.9148; 0.02; 0.0003it's only a coefficent 's change of the center's equation in degrees(we just have to multiply the first coefficients by 180/π). wiki explains how to get the center's equation.

280.4665= it's the ecliptic's longitude value the 01/01/2000 to 12h.

0.98564736= it's the average angular's velocity about earth around the sun on a tropic year, So we find this value by calculating 360/365.2421898.

-2.46569; 0.053; 0.0014= as the center's equation it's a switch from the first coefficients to coefficients in degrees.

b) declination angle

As we have seen in “programming a TI89 for navigation”, there is a lot of several ways to get sun declination; by using methods with less or more parameters, thus less or more accuracy and less or more difficulties to program it too… not easy. So we will remove "useless" data for our case. we will not take care about planets and moon leverage in our calculations, but Kepler's laws used on the earth/sun system only, it will be not bad at all to begin.

the easiest formula to have declination is the formula which have a simple periodic function's framework (sinus, cosinus) like Desmond Fletcher's equation:

D=23.45*sin(2π/365(284+j)) avec j jours juliens. Here 23.45 answers to the sinus function range, thus to the highest value reachable by declination (during solstices). 2π/365 matches to the sinus function's period, here it will repeat itself every 365 days. 284 is a function's shift. this formula doesn't have in charge the earth ellipticity's trajectory and consider its motion as a cercle's shape. the easiest method but really inaccurate.

for declination an another way to do would be to write a periodic function but by parts, each part matches for a season, thus a function with 4 parts (see “programming a TI89 for navigation”). it works well, easy to do because all we have to know is solstices and equinoxes dates for the following year, however its weakness is that it works just for one year. in fact from a year to an another equinoxes and solstices hours change by following an average value, but the truth is that there are gaps with ranges of few minutes, and at the end this differences cumulation brings…highest differences!

on internet it's easy to find formulas for declination given by IMCCE, efficient for a long or short period. we could also deduce it from the time's equation because we know that:

sin(declination)= sin(ecliptic's obliquity)*sin(sun's true longitude).

Why ? because , to prove it it's almost like method used for the rectangular and spheric triangle (allow us to find formulas for azimut and calculated height for a star, planet...), as i am lazy and it's very well explained on the following website page : https://sites.google.com/site/astronomievouteceleste/4—mouvement-du-soleil

from all this we could deduce declination's formula to write after our program's beginning:

So for declination, 3 programmation lines are enough when we know AHG. before going ahead we will check efficiency thanks to the IMCCE website, just to notice if it's accurate enough.

Test

by taking 10 several values ​​we get that:

the 15/12/2019 at 2:22 p.m. : -23°16’13.91”/36°44’27.88”→-23°16'/17h31

the 2/11/2043 2h : -14°40’9.90”/214°6’50.43”→-14°40'/14h29

the 28/2/2055 at 4:47 p.m.: -7°47’20”/68°38’28.54”→-7°47'/22h46

the 12/5/2128 at 5h06: 18°28’20.7”/257°22’25.4”→18°13'/3h18

the 7/07/2184 at 21h12: 22°17’33.11”/136°34’51.76”→22°23'/7h12

the 18/10/2079 23h10: -9°56’58.8”/171°15’1.75”→-9°57/13h35

the 13/6/2046 at 12h11: 23°14’12.71”/2°42’14.74”→23°14'/5h28

the 6/09/2033 at 18h55: 6°6’8.13”/104°12’36.45”→6°6'/11h03

the 1/03/2085 at 4h22: -7°19’47.17”/242°27’31.56”→-7°20'/22h51

the 19/02/2074 at 3h18: -11°10’21.23”/226°4’8.9”→-11°10'/22h12

on the left it's our prog results for declination and greenwich hour angle, and on the right it's the IMCCE website's results “calcul des éphémérides” they give them in declination and right ascension. not to bad it seems, we begin to lose accuracy after 2080 on declination so we have time… we know that right ascension is the angle between vernal point's meridian and a planet's one, measured from east and written in hours, and it's the 360 degrees complement of verse ascension. but we have already found a formula able to calculate vernal point's meridian in our prog to get our position with stars :

100+5.5/60+(l-2451545-5843.5)*(360/365.2421898+360)-0.0524*16/(24*3600*36525)-0.0524*(l-2451545-5843.5)/(24*3600*36525)=vernal point's angle. thus by putting it after in our prog it gets that:

when we compare it with right ascensions given by IMCCE it's not too bad. we could say that until 2080 roughly it's accurate, and after it's decreasing in accuracy. we could consider that our ephemeris part is OK for navigation. we will be able to challenge the next step that's cool. for this part we could have been more meticulous by integrate some stable parameters used like variables like excentricity, duration of tropic and anomalistic year, and inclination angle's earth on ecliptic. however for our needs it will be enough i guess.

3) intercept and true azimut calculations

now we have declination and AHL calculated at watch time, the next level will be to define the calculated height (trigo formula like stars), true height (it's the instrumental height with corrections), and thus intercept and the true sun's azimut.

corrections to add up with instrumental height:

– Refraction of light (depends on the height above the horizon line).

– parallax (depends on the radius of earth, earth/sun distance and height above the horizon's line).

– sun's half diameter (roughly 16 minutes).

– eye's watcher height above the sea's level.

– collimation.

our program draft gives this:

this program allows to calculate meridians too, you just have to push 1 or 2 if you are back to the south or north when you are focusing the Sun.

Attenzione: test !

i have checked the device with navastro's webpages (“navastro.free.fr/meridienne.htm” to do a meridian and “navastro.free.fr/exos/exercices.htm” for straight heights). results are not totally the same but it remains "not bullshit" enough. A little disappointed, I dig it and notices an error in the program above; at the AHL calculation line (ahl) instead of having:

: mod(qq-rrr, 360)→ ahl in fact it's :mod(qq+rrr,360)→ahl

And for intercept instead of having:

: hv-hc→in in fact it's (hv-hc)*60→in

The final program that works is that:

By making these changes, this time we find almost the same results compared with navastro! well, lets Haddock's captain concludes with philosophy…